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Abstract 

A  model  of  a  polymer  electrolyte  fuel  cell  (PEFC)  is  developed  that  captures  dynamic  behaviour  for  control  purposes.  The  model  is  mathematically 
simple,  but  accounts  for  the  essential  phenomena  that  define  PEFC  performance.  In  particular,  performance  depends  principally  on  humidity, 
temperature  and  gas  pressure  in  the  fuel  cell  system.  To  simulate  accurately  PEFC  operation,  the  effects  of  water  transport,  hydration  in  the 
membrane,  temperature,  and  mass  transport  in  the  fuel  cells  system  are  simultaneously  coupled  in  the  model.  The  PEFC  model  address  three 
physically  distinctive  fuel  cell  components,  namely,  the  anode  channel,  the  cathode  channel,  and  the  membrane  electrode  assembly  (MEA).  The 
laws  of  mass  and  energy  conservation  are  applied  to  describe  each  physical  component  as  a  control  volume.  In  addition,  the  MEA  model  includes 
a  steady- state  electrochemical  model,  which  consists  of  membrane  hydration  and  the  stack  voltage  models. 

©  2007  Published  by  Elsevier  B.V. 
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1.  Introduction 

Fuel  cells  may  be  used  to  replace  the  conventional  inter¬ 
nal  combustion  engine  in  vehicles  because  of  their  capacity  to 
achieve  a  high  efficiency.  In  particular,  they  are  being  considered 
as  suitable  potential  power  sources  for  zero-emission  vehicles 
(ZEVs).  Most  automobile  manufacturers  consider  the  fuel  cell 
electric  vehicle  to  be  logical  successor  to  the  hybrid  electric  vehi¬ 
cle,  and  major  automotive  companies  have  already  presented 
cars  that  use  fuel  cells  as  the  power  plant.  Among  the  various 
types  of  fuel  cells,  the  polymer  electrolyte  fuel  cell  (PEFC)  has 
the  greatest  potential  for  use  in  automotive  applications  because 
of  its  low  operating  temperature,  solid-type  electrolyte,  high 
efficiency,  and  high  power  density  [1].  PEFC  technology  has 
now  reached  the  testing  and  demonstration  phase,  and  engineer¬ 
ing  development  and  optimization  are  vital  for  realization  of 
commercial  products  [2].  It  is  especially  necessary  for  PEFC 
technologies  to  be  developed  in  tandem  with  dedicated  control 
systems  which,  in  turn,  requires  effective  modelling  of  technol¬ 
ogy.  The  operation  and  physical  phenomena  within  the  PEFC 
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system  can  be  clearly  understood  through  this  modelling  pro¬ 
cess.  The  knowledge  acquired  through  the  modelling  the  PEFC 
and  simulations  using  this  model  can  also  be  applied  to  the 
design  of  control  systems  for  the  PEFC. 

2.  Model  setup  and  assumption 

Although  the  PEFC  model  developed  in  this  study  is 
mathematically  simple,  it  accounts  for  the  essential  phenomena 
that  define  PEFC  performance.  The  performance  of  the  PEFC 
depends  principally  on  the  properties  of  the  membrane  and 
the  thermodynamics  of  the  gas  in  the  fuel  cell  system,  and 
the  model  emphasizes  these  elements  accordingly.  The  PEFC 
system  consists  of  the  three  physically  distinct  parts  shown  in 
Fig.  1,  namely,  the  anode  channel,  the  cathode  channel,  and 
the  membrane  electrode  assembly  (MEA)  [3].  The  mass  and 
energy  conservation  laws  are  applied  to  describe  each  physical 
component  as  a  control  volume.  In  addition,  the  MEA  model 
contains  a  steady- state  electrochemical  model  that  consists  of  a 
membrane  hydration  model  and  a  stack  voltage  model  [4] .  This 
model  is  the  lumped  parameter  model,  so  spatial  variations 
of  the  temperature,  pressure,  and  concentration  are  ignored. 
Dynamics  for  flow  and  temperature,  however,  are  included.  The 
inputs  are  electric  current,  mass  flow  rate,  humidity,  and  the 
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Nomenclature 

a 

activity  of  water 

A 

fuel  cell  active  area  (cm2) 

can  ,v 

anode  water  concentration 

Cc  A,V 

cathode  water  concentration 

CP 

specific  heat  at  constant  pressure  (kJkg-1  K) 

Cy 

specific  heat  at  constant  volume  (kJkg-1  K) 

£>w 

diffusion  coefficient  of  water  concentration 

£CV 

internal  energy  of  control  volume  (kJ) 

^■Nerst 

reversible  open  circuit  voltage  (V) 

F 

Faraday  constant  (C) 

A/feact 

enthalpy  of  reaction  (kJkmol-1) 

(M) 

heat  convection  coefficient  (kJ  K-1) 

1 

stack  current  (A) 

i 

stack  current  density  (A  cm-2) 

k 

effective  permeability  (cm2) 

^sat 

permeability  (cm2) 

K 

nozzle  coefficient  (kg  s-1  atm) 

m 

mass  (kg) 

m 

mass  flow  rate  (kg  s-1) 

M 

molecular  mass  (kgkmol-1) 

N 

molar  flux  (kmol  m_2s“) 

^cell 

number  of  cell 

nd 

electro-osmotic  drag  coefficient 

p 

pressure  (kPa) 

Q 

rate  of  heat  transfer  (kW) 

T 

temperature  (K) 

To 

standard  temperature  (K) 

tM 

thickness  of  the  membrane  (cm) 

V 

channel  volume  (cm3) 

V 

molar  volume  (cm3) 

w 

electric  power  (kW) 

X 

mole  fraction 

Greek  letters 

7?  act 

activation  loss  (V) 

0  ohmic 

Ohmic  loss  (V) 

Oconc 

concentration  loss  (V) 

X 

water  content 

PM, dry 

membrane  dry  density  (kg  cm-3) 

&  m 

membrane  conductivity 

Subscripts 

AMB 

ambient  condition 

AN 

anode 

CA 

cathode 

cone 

concentration 

d 

electro-osmotic  drag 

h2 

hydrogen  gas 

L 

liquid  type  of  water 

m 

membrane 

M 

membrane  electrode  assembly 

n2 

nitrogen  gas 

o2 

oxygen  gas 

ohmic 

Ohmic 

prod 

production 

sat 

saturation  condition 

V 

vapor  type  of  water 

W 

water 

Superscripts 

ave 

average 

conv 

convection 

diff 

diffusion 

I 

inlet 

latent 

latent 

0 

outlet 

osmotic 

electro-osmotic  drag 

react 

reaction 

trans 

translation 

temperatures  of  the  supplied  gas  and  coolant.  Some  states  which 
cannot  be  directly  calculated  by  input  variables,  such  as  the 
pressure,  humidity,  concentration  and  temperature  of  gas  in  the 
channel,  are  determined  by  the  dynamics  of  the  fuel  cell  system. 
The  output  variables  of  a  fuel  cell  system  such  as  exhaust 
temperature,  exhaust  mass  flow  rate,  and  the  stack  voltage  are 
considered  to  be  functions  of  the  states  and  the  input  variables. 
The  main  assumptions  of  this  PEFC  model  are  as  follows: 

•  the  gases  are  ideal  and  well-mixed  in  the  channel, 

•  there  is  no  pressure  gradient  in  the  channel, 

•  there  is  no  temperature  gradient  in  the  channel, 

•  there  is  no  current  gradient  across  the  membrane, 

•  the  MEA  temperature  is  uniform, 

•  the  mole  fractions  in  the  channel  and  exhaust  gas  are  equal, 

•  the  temperatures  of  the  condensed  water  and  the  MEA  are 
equal, 

•  the  radiation  heat- transfer  effect  is  negligible. 

3.  Mathematical  model 

3. 1.  Review  of  thermodynamics 

3.1.1.  Saturation  vapour  properties 

Water  is  always  accompanied  with  the  dry  gas  in  the  chan¬ 
nel  because  of  hydration  of  the  membrane  and  the  generation 
of  water  due  to  the  chemical  reaction.  So  the  phase  of  water 
in  the  PEFC  is  important  factor  in  determining  the  thermody¬ 
namic  states  such  as  pressure,  and  temperature.  The  saturation 
vapor  pressure  pSat(T)  (kPa)  according  to  vapour  temperature  is 
modelled  using  the  table  of  saturated  water  [5]. 

logio Psat(T)  =  -20.92  +  0.1437  -  3.39  x  10~4T2 
+3.85  x  l<r773  -  1.69  x  10“1074 


(1) 
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Fig.  1.  PEFC  model  structure. 


where  T  is  the  vapour  temperature  (K).  The  mass  of  saturation 
vapour  is  determined  by  the  Ideal  Gas  Law. 


my,sat 


Psa\V 


M^2oRT 
The  derivative  of  mysat  is  given  by: 

dmv?sat 


d  t 


3mysat  dT  3m  v;  sat  dPsat 

3  T  d  t  3/?sat  d  t 


(2) 


(3) 


where  3mysat/3r  =  — ( pSSLty/R)T~ 2  and  3mysat/3pSat  =  V/(RT). 
The  above  relationship  requires  the  determination  of  the  time 
derivative  of  psat(T). 


d/?sat 

d t 


Psat  In  10 


dr  d(log10psat) 


d  t 


dr 


(4) 


3.2.1.  Anode  channel  flow  model 

The  mass  continuity  of  hydrogen  and  water  is  used  to  develop 
the  flow  model  in  the  anode  control  volume. 


dmAN,H2  .  i 

-  "‘an, *2 


■  m 


•  O 


an,h2 


m 


react 

H2 


dm  an, 


w 


d  t 


=  m 


•  i 


AN,V  ' 


■  m 


•  O 


(6) 


AN,V 


■  m 


•  trans 


The  water  inside  the  anode  control  volume  is  one  of  phases, 
vapour  or  liquid,  depending  on  the  amount  of  maximum  vapour 
mAN,v,sat  [5],  which  is  in  turn  determined  by  the  saturation  pres¬ 
sure  pan, v, sat  corresponding  to  the  temperature  of  the  anode 
gas. 

If  mAN,W  >  mAN,V,sat> 


3.1.2.  Latent  heat 

The  latent  heat  effect  is  a  crucial  factor  that  affects  the  heat 
transfer  in  the  PEFC  system,  because  a  large  quantity  of  water 
is  generated  by  chemical  reactions  and  is  changed  between  the 
liquid  and  vapour  phases  according  to  the  dew  point  and  the 
vapour  pressure.  The  heat  of  evaporation  or  condensation  //vapor 
may  be  estimated  as  a  function  of  temperature  by  following  Eq. 
[6]: 

Hv apor  =  45,  070  -  41.97’  +  3.44  x  10“37’2 

+2.54  x  10“67’3  -  8.98  x  10_107’4  (5) 


3.2.  Anode  channel  model 

Hydrogen  and  vapour  enter  the  anode  channel  (Fig.  2),  some 
of  which  diffuses  through  the  gas-diffusion  layer  (GDL)  and 
some  of  which  leaves  the  channel  [4].  Heat  transfer  from  the 
MEA  to  the  gas  in  the  channel  affects  the  temperature  and 
pressure  of  the  gas  in  the  anode  channel.  The  anode  chan¬ 
nel  is  considered  as  a  control  volume  in  describing  these 
phenomena. 


mAN,V  =  mAN,V,sat 

^AN,V  =  mAN,V,sat 

mAN,L  =  mAN,W  mAN,V 
.  vapor 

mAN,L  =  ~mAN  v  =  mAN,w  ^an,v 


(7) 


mAN.H,  +,^AN.V 


Fig.  2.  Anode  channel  model. 
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Or  if  mAN,w  <  mAN,v,sab 


^AN,V  =  ^AN,W 
^AN,V  =  ^AN,W 

r\  (°) 

mAN,L  =  0 

.  vapor  ~ 

PlAN,L  —  — ^AN,V  —  ^ 

The  mass  flow  rate  of  the  anode  exhaust  gas  m^N  is  assumed 
to  be  proportional  to  the  pressure  difference,  i.e., 

^an  =  kan(pan  7>amb)  (9) 


It  is  assumed  that  the  mole  fractions  of  the  outlet  flows  are 
the  same  as  those  of  the  gas  in  the  anode  control  volume.  Using 
anode  channel  pressure  and  mass  fractions,  hydrogen  flow  and 
water  flow  at  the  outlet  can  be  given  by: 


.  o  _  MiXANjKANiPAN  ~  7>AMb)  . 


mAN,z  — 


J2iMjxAN,j 


i  =  H2,  H20 


(10) 


Faraday’s  law  is  used  to  calculate  the  mass  flow  rate  of 
hydrogen  consumption  in  reaction  [7].  The  mass  flow  rate  is 
determined  according  to: 


■  react 
%2 


/ 

=  Mu2nc& n  — 


(ID 


3.2.2.  Anode  thermal  model 

The  energy  rate  balance  for  the  anode  control  volume  is: 


dE' 


cv 

AN 


d  t 


—  ^AN 


+  ^an^an  -  ,h^K]h 


O  ho 

AN'7  AN 


(12) 


where  Q  anV  is  the  heat  transfer  rate  across  the  MEA  to  the  anode 

channel,  and  ^an^an  —  ^an^an  *s  i^e  net  enerSY  rate 
by  the  mass  flow.  The  convective  heat  transfer  Q™ is  caused 
by  the  temperature  difference  between  the  MEA  and  the  anode 
[3,8,9]: 

Ganv  =  (^A)an(Im  -  Tan)  (13) 


where  (M)an  is  the  heat  convection  coefficient  of  the  anode 
channel.  The  net  energy  rate  by  the  mass  flow  is  given  by: 

^  ^AN^AN  ^AN^AN 

=  ^an,h2cp,H2(Tan  —  To)  —  ^an,h2cp?h2(Tan  —  To) 

— 777  an,h2  cp,  h2  (  Tan  —  To)  +  ^an,vcp,v(Tan  —  To) 

mAN,VcP,v(TAN  To)  +  mAN,Vcvy^  AN,V  To) 

— myanscP;v(TA^ny  -  Tb)  (14) 

where  Tj^ y  represents  the  temperature  of  the  vapourized  water 
or  the  temperature  of  the  condensed  water  in  the  anode  channel, 
and  T^ny  represents  the  temperature  of  the  vapor  moving  from 
the  MEA  to  the  anode  channel  or  from  the  anode  channel  to  the 
MEA.  r;aNP°vr  and  are  determined  by  the  sign  of  m^y 

and  myans.  If  m^y  is  positive,  Tj^ y  is  the  same;  otherwise 
An,v  is  Tan-  An’v  is  TAN  if  m%ms  is  positive  and  TM  if  mlyim 


is  negative.  The  internal  energy  of  the  anode  control  volume  can 
be  expressed  as  follows: 

=  (^an,h2G,h2  +  ^an,vCv,v)(Tan  To)  (15) 

Using  Eqs.  (12)  and  (15),  the  time  rate  of  the  temperature  for 
anode  control  volume  can  be  expressed  as  follows: 

dTAN  _  -&AN  (^AN,H2^v,H2  +  ^AN,V^v,y)(TaN  Tq) 
d  t  ^an,h2G,h2  +  ^an,vG,v 

(16) 


3.3.  Cathode  channel  model 

The  mass  flow,  energy  flow,  and  gas  pressure  in  the  cathode 
control  volume  are  similar  to  those  in  the  anode  control  volume. 
Flows  of  oxygen,  nitrogen  and  water  in  the  cathode  are  consid¬ 
ered  as  in  Fig.  3.  Oxygen  molecules  are  involved  in  the  chemical 
reaction  and  move  to  the  membrane  and  react  with  protons.  As 
a  result,  liquid  water  is  produced  at  the  catalyst  of  the  cathode 
side  of  the  MEA  [10,11].  Moreover,  heat  transfer  by  temperature 
gradient  to  the  MEA  affects  the  temperature  and  pressure  of  the 
gas  in  the  cathode  channel. 

3.3.1.  Cathode  flow  model 

The  cathode  flow  model  describes  the  flow  of  air  and  vapour 
inside  the  cathode  channel.  The  mass  conservation  law  and  the 
Ideal  Gas  Law  are  used  to  determine  the  pressure  of  the  gas  in 
the  cathode  control  volume. 

'“CA.Cb  =  '«CA,02  -  '«CA,02  -  ™0 T 

«CA,N2  =  "*CA,N2  —  "*CA,N2  (17) 

.  .  t  •  O  -  trans  1  •  prod 

^CA,W  =  ^CA,V  mCA,V  +  mV  +  mL 

The  behaviour  of  the  water  inside  the  cathode  channel  is 
somewhat  different  from  that  or  of  the  water  inside  the  anode 
channel  due  to  the  production  of  water  reacting  at  the  cathode 
side  of  the  membrane  close  to  the  GDL.  It  is  assumed  that  all 
of  the  water  produced  by  the  chemical  reaction  is  liquid.  The 
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water  inside  the  cathode  channel  is  in  one  of  two  phases,  i.e., 
vapour  and  liquid,  and  depends  on  mcA,v,sat?  which  is  in  turn 
determined  by  the  saturation  pressure  [5]. 

If  mcA,W  >  mCA,V,sat> 


PICK  ,V  =  ^CA,V,sat 
PlCAy  =  ^CA,V,sat 

^ca,l  =  mcA,w  mcA,v  (18) 

PI  CA.L  =  ^CA,W  -  PlCAy 

.  vapor  .  .  .  prod 

PI  CA .  V  —  mCA,V  +  ^CA.L  mCA,W 

Else  if  mcA,w  <  mcA,v,sat, 


pic  Ay 
pic  Ay 
PI  CA,L 
^CA,L 


vapor 

^CA,V 


=  PI  CA,W 
=  }P  CA .  W 
=  0 


=  PIC  A,W  PlCAy 

.  prod 
=  mCA,L 


0 


(19) 


It  is  assumed  that  the  mole  fractions  of  the  outlet  flows  are  the 
same  as  those  of  gas  in  the  cathode  control  volume.  The  cath¬ 
ode  outlet  flow  of  each  species  can  be  expressed  using  pressure 
difference  and  mass  fractions  [3]: 


m 


.  o  _  KcAMiXcA,i(pCA  —  Pamb)  . 


CA,/ 


72iMix  CA,; 


i  =  02.  N2,  H20  (20) 


The  mass  flow  rate  of  oxygen  for  the  chemical  reaction  is 
determined  by  the  stack  current  [7],  i.e., 


mgf  =  Mo^cell  — 


(21) 


3.3.2.  Cathode  thermal  model 

The  energy  rate  balance  of  the  cathode  control  volume  is 
given  by: 


d  t 


QbT  +  y^cAh 


I 

CA 


L«ca^ca 


(22) 


where  Qqa  ls  the  rate  of  heat  transfer  across  the  MEA  to  the 

cathode  channel,  ^ca^ca  —  ^ca^ca  is  l^e  net  enerSY 
rate  by  mass  flows.  Heat  convection  QqA  Is  caused  by  a  dif¬ 
ference  in  the  temperature  between  the  MEA  and  the  cathode 
[3,8,9]: 


Qca  =  (hA)CA(TM  -  TCa) 


(23) 


Fig.  4.  Membrane  electrode  assembly  model. 


_^ca,n2cp,n2(Tca  -  To)  +  ^ca,vcp,v(Tca  _  To) 
-mgA,vCp,v(7cA  -  T0)  +  mfnscpy(T^  -  T0) 
+-cTv^v(^caPv  -  To)  (24) 


where  T™A  y  represents  the  temperature  of  the  vapourized  water 
or  that  of  the  condensed  water  in  the  cathode  channel,  and  TqA- y 
represents  the  temperature  of  the  vapour  moving  between  the 
MEA  and  the  cathode  channel  or  from  the  cathode  channel  to 
the  MEA.  T^°y  and  T^A  y  is  determined  by  the  signs  of  m™A  y 
and  myans.  If  m™A  y  is  positive,  T™A  y  is  same  as  Tm •  Otherwise, 
TqA  y  is  7ca-  T^ny  is  7m  if  myans  is  positive  and  Tqa  if  Wyans 
is  negative.  The  internal  energy  of  the  cathode  control  volume 
can  be  expressed  as  follows: 

cv 

ECA  —  {(mCA,02Cv,02  +^CA,N2Cy,N2  +  PiCAycvy) 

x(Tca-To)}  (25) 


Using  Eqs.  (22)  and  (25),  the  time  rate  of  the  temperature  for 
cathode  control  volume  can  be  expressed  as  follows: 


dTcA  _  TgX  -  (mcA,02Cy,Q2  +  mCA,N2Cy,N2  +  PiCAycYy)(TcA  ~  Tp) 

d  t  pica,o2cv,o2  +  ^ca,n2cv,n2  +  PiCAy  cYy 


(26) 


where  (M) ca  is  the  heat  convection  coefficient  of  the  cathode 
channel.  The  net  energy  rate  by  mass  flows  is  given  by: 

53  ^CA^CA  -  53  ^CA^CA 

=  "1ca,o2cp,o2(7cA  —  To)  —  ^ca,o2cp,o2(7ca  —  To) 
-'«CACt02cP,02(7cA  -  To)  +  ™CA,N2cP,N2<Tca  -  7o) 


3.4.  Membrane  electrode  assembly  model 

Three  sub-models  included  in  the  MEA  model  are  the  mem¬ 
brane  hydration  model,  the  stack  voltage  model,  and  the  MEA 
thermal  model.  Simplified  features  of  the  MEA  are  shown  in 
Fig.  4,  which  illustrates  the  energy  and  the  mass  balances  within 
the  MEA  control  volume. 
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3.4.1.  Membrane  hydration  model 

Water  molecules  are  dragged  across  the  membrane  from  the 
anode  to  the  cathode  by  the  protons.  This  phenomenon  is  called 
electro-osmotic  drag.  There  is  a  water  concentration  gradient 
across  the  membrane  caused  by  the  difference  in  the  humidity 
between  the  anode  and  cathode  flows  and  which  in  turn  causes 
diffusion.  The  pressure  gradient  across  the  membrane  causes  the 
convection  within  the  pores  of  the  membrane.  By  combining  the 
three  water  transports,  the  water  flow  across  the  membrane  m yans 
can  be  written  as  follows  [5,10-14]: 


•  trans _ 

riy  — 


(TCA,V  —  cAN,v) 


m^ans  =  ft cen  X  A  X  Mh20  X  ftyans 


k  (Pc a  ~  Pan) 
ll  Vo  tM 


(27) 


where  A  is  the  active  area  of  the  fuel  cell.  The  membrane  water 
content  and  the  electro-osmotic  and  diffusion  coefficients  can 
be  calculated  using  the  activities  of  the  water  in  the  anode  and 
the  cathode  [4,5].  It  is  assumed  that  the  temperature  of  water 
adjacent  to  the  membrane  is  the  same  as  that  of  the  MEA.  So 
the  activities  of  water  at  the  anode  and  the  cathode  sides  of  the 
membrane  are  determined  by: 


_  *an,wPan  _  _  -^ca,w  Pc  A 
°AN  ~~  />sat(7 m)  ’  aCA  ~  Psat(7 m) 

_  CL  AN  +  CLC  A 
^ave  —  Z 


(28) 


permeability  and  the  water  volume  fraction  in  the  membrane. 
The  effective  permeability  can  be  expressed  by  [13]: 


k  —  ^sat 


(32) 


where  /l  is  the  maximum  water  volume  fraction  and  /  is  the 
water  volume  fraction.  The  water  volume  fraction  is  expressed 
by  [13]: 


A  Vw 

Vm  +  ^w’ 


(33) 


where  Vm  is  the  partial  molar  volume  of  the  membrane  and  Vw 
is  the  molar  volume  of  water  [13]. 


Vw  = 


Mw 

Pw 


37m,  dry 
PM,  dry 


(34) 


3.4.2.  Stack  voltage  model 

The  voltage  is  a  function  of  stack  current,  cathode  pres¬ 
sure,  reactant  partial  pressure,  MEA  temperature,  and  membrane 
humidity.  The  cell  voltage  is  determined  by  the  combination  of 
T^Nerst,  ?7act,  77ohmic,  T] COnc  [5,15—16]: 

l^cell  =  ^Nernst  —  Pact  ~  ^?ohmic  —  Oconc  (35) 

The  voltage  ^Nernst  is  the  reversible  open-circuit  voltage,  which 
is  derived  from  the  change  in  Gibbs  free  energy  [4,5,7]. 


The  membrane  water  content  A/  is  calculated  from  water 
activities  ft/  [5]: 


*i  = 


0.043  +  17.81ft/  —  39.85ft?  +  36.0ft?,  0  <  ft/  <  1 

14  +  1.4(ft/  —  1),  1  <  ft/ 


^Nemst  =  1-229  -  (8.5  X  10"4)(7m  -  298.15) 

+(4.308  x  10_5)7M(ln /?an,h2  +  0.5  ln/?cA,o2) 


(36) 


i  =  AN,  CA 


(29) 


The  average  membrane  water  content  Aave  is  the  arithmetic 
mean  of  the  water  content.  The  electro-osmotic  drag  coefficient 
ftd  and  the  water  diffusion  coefficient  Z>w  are  calculated  from 
Aave  as  follows  [4,5]: 

2.5  f 

ftd  =  7^- Aave,  E>w  =  P>x  exp  (  2416 

(30) 


The  water  concentrations  at  the  membrane  surfaces  on  the 
anode  and  cathode  sides  are  functions  of  the  membrane  water 
content  [5]: 


CAN,V 


PM,  dry 
37m,  dry 


A  AN, 


CCA,V 


PM,  dry 

37m,  dry 


Aca 


(31) 


For  the  liquid  water  flow  through  the  membrane,  Darcy’s  law 
is  used  to  express  an  influence  of  the  pressure  gradient.  In  Eq. 
(27),  the  water  flux  caused  by  the  pressure  gradient  is  propor¬ 
tional  to  the  effective  permeability  and  the  inverse  of  viscosity 
of  water.  The  effective  permeability  is  a  function  of  the  absolute 


Activation  losses  are  caused  by  the  need  to  move  electrons 
and  to  break  and  form  chemical  bonds  at  the  anode  and  cathode 
[5].  The  activation  losses  are  described  by: 

Oact  =  vo  +  va{\  ~  e~c^)  (37) 

Ohmic  losses  are  caused  by  the  resistance  of  the  membrane 
to  the  transfer  of  protons  and  resistance  to  the  flow  of  electrons 
through  the  electrode  material  and  various  interconnections.  The 
resistance  due  to  the  movement  of  electrons  is  ignored  in  this 
model.  The  Ohmic  losses  are  proportional  to  the  current  density 
[5]: 


^ ohmic  =  —U  =  (b\\ Am-fel2)  exp  (  b2 

(38) 

where  crm  is  the  membrane  conductivity.  Concentration  losses 
result  from  the  reduction  in  concentration  caused  by  insufficient 
quantities  of  reactant  to  the  electrode  surface  [5],  i.e., 

Oconc  =  i(c2  T~  ^  ,  *max  =  2.2  (39) 

\  *max  / 


1  1 
303  _  7k 
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3.4.3.  ME  A  thermal  model 

The  energy  rate  balance  of  the  MEA  control  volume  can  be 
expressed  as: 

A  rCV 

^  =  Qm  +  ^mCC  +  “  22  (4°) 

where  Qy^  is  the  rate  of  heat  transfer  within  the  MEA,  W^c  is 
the  electric  power  generated  by  the  PEFC,  is  the  thermal 

energy  transfer  rate  due  to  flows  coming  into  the  MEA,  and 
1S  thermal  energy  transfer  rate  due  to  flows  going  out 
of  the  MEA.  The  heat  transfer  rates  in  the  MEA  consist  of  the 
seven  terms  [3,6,8,11]: 

fiCV  _  Areact  /Sconv  Aconv  Aconv  Ajconv 
—  ^AN  —  ^CA  —  *3 CL  ~  ^AMB 

-  Gan"'  -  Gc!T  (41) 


Ai react  -react  A  ri-react  //io\ 

Qm  =  mAN,U2AH  (42) 

where  AEEQiiCt  is  the  lower  heating  value  of  the  fuel  [3,8-9].  The 
convective  heat  transfer  to  ambient  air  and  to  coolant  water  is: 


Game  =  (hA)AMB(TM  -  Tame)  (43) 

GcT  =  (hA)CLATlm  =  mcLCp.L^CL  -  T]CL)  (44) 


where  (M)cl  and  (hA) amb  are  the  heat  convection  coefficients 
of  the  cooling  water,  A  Tim  is  the  log  mean  temperature  differ¬ 
ence  between  coolant  inlet  and  coolant  outlet,  and  m<x  is  the 
mass  flow  rate  of  the  coolant  [3,8-9].  The  log  mean  temperature 
difference  between  coolant  inlet  and  outlet  is  defined  as  follows: 


A  Tim  = 


A t£l  -  AT, 


CL 


ln(AT°L/AT)L) 


(45) 


ATS.  =  Tm-T, 


CL’ 


A7£l  =  Tm-T( 


CL 


(46) 


When  liquid  water  exists  in  the  channel,  it  tends  to  be  on 
the  channel  surface.  Since  the  heat  transfer  rate  between  the 
liquid  water  droplets  and  the  MEA  is  always  greater  than  that 
between  the  liquid  water  droplets  and  the  gas,  it  is  assumed 
that  the  temperature  of  the  liquid  water  and  the  MEA  are  equal. 
Therefore,  it  is  assumed  that  the  phase  change  occurs  on  the 
surface  of  the  channel  [6].  The  latent  heat  transfer  due  to  the 
phase  change  of  water  is  given  by: 

.vapor  rr  vapor  .vapor  ry  vapor 

Ailatent  _  mAN  n  /^latent  _  mCA  n  <ah\ 

CaN  -  Mh,o  '  QcA  -  M„,  O  (47) 

The  term  Wyfc  represents  the  electric  power  converted  from 
the  electrochemical  reaction  with  nce\\  single  cells,  each  of  area 
A  [3]. 

K\ec  =  -  Kell  X  Been  x  /  X  A  (48) 


The  energy  transfer  due  to  mass  transfer  across  the  boundary 
of  the  MEA  is  expressed  as: 

=  mrI|2actCp,H2(7AN  To)  +  mr0e2actCp,°2(TcA  To) 

+m yans cp? v ( ^an!v —To)—  mfn*cpy(T^  -  T0)  (49) 


Table  1 

Thermodynamic  constants 


Parameter 

Value 

Parameter 

Value 

R 

8.314  kJkmol_1K 

cv,02 

0.668  kJ  kg”1  K 

F 

96,485 

cv,N2 

0.744  kJ  kg”1  K 

cp,h2 

14.472  kJ  kg-1  K 

Cv,V 

1.419kJkg->  K 

Cp,02 

0.928  kJ  kg- XK 

Mh2 

2.016  kg  kmol-1 

cp,n2 

1.041  kJ  kg"1  K 

Mo2 

32kgkmol_1 

cp,v 

1.881  kJ  kg"1  K 

MN2 

28  kg  kmol-1 

cp,L 

4.184  kJ  kg"1  K 

4^h2o 

18.02  kg  kmol-1 

cv,H2 

10.302  kJ  kg"1  K 

A  LTeact 

285.84  kJ  mor1 

The  time  rate  of  change  of  the  MEA  temperature  is: 

_  6mV  +  ^mCC  +  X] 

d  t  pimcv,m 

Parameters  adopted  for  thermodynamic  properties  of  species 
are  summarized  in  Table  1 .  The  specific  heats  of  species  and  gas 
constants  are  collected  from  some  published  papers. 

4.  Simulation  results  and  discussion 

The  PEFC  stack  is  a  multi-input  and  -output  system.  The 
inputs  of  the  model  are  mass  flow  rate,  relative  humidity,  and 
inlet  gas  temperature.  The  flow  rate  and  temperature  of  the  cool¬ 
ing  water  are  the  inputs  of  the  model  that  determine  the  operating 
conditions.  Simulation  results  show  the  effects  of  varying  the 
inputs  of  the  stack  system.  In  this  section,  the  simultaneous 
effects  of  temperature,  humidity  and  pressure  are  elucidated 
using  simulation  results.  Every  operating  condition  except  for 
electric  load  current  is  fixed.  The  relative  humidity  of  each 
channel  is  set  to  100%  and  the  flow  rate  and  temperature  of 
cooling  water  are  set  to  50  L  per  min  and  40  °C,  respectively. 
The  hydrogen  gas  feed  is  fixed  to  1.371  x  10-4gs-1  which 
is  1.5  times  the  stoichiometric  hydrogen  flow  rate  correspond¬ 
ing  to  a  current  load  of  250  A.  The  air  feed  for  the  cathode 
is  set  to  6.539  x  10-3gs-1,  i.e.,  twice  the  stoichiometric  air 
flow  rate  corresponding  to  a  current  load  of  250  A.  The  anode 
inlet  and  cathode  inlet  gases  are  set  to  65  °C  and  are  humidified 
to  a  relative  humidity  of  100%  corresponding  to  the  inlet  gas 
temperature. 

Table  2  shows  the  parameters  used  for  the  simulation.  These 
parameters,  such  as  cell  count,  channel  volume,  and  heat  con- 


Table  2 

Parameters  dependent  on  the  stack  design 


Parameter 

Value 

References 

TTcell 

35 

[8] 

tm  (cm) 

0.0175 

[4] 

Van  (m3) 

0.005 

[5] 

Vca  (m3) 

0.01 

[5] 

(hA) an  (kWK-1) 

2  x  10"3 

[8] 

(hA) ca  (kWr1) 

10  x 10-3 

[8] 

(M)amb  (kWK"1) 

17  x 10“3 

[8] 

(hA) CL  (kWr1) 

181.8  x  10-3 

[8] 

mMcp,M  (kJK-1) 

35 

[3] 
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Fig.  5.  Electric  load  current. 

vection  coefficients  are  related  to  stack  design,  and  they  may  be 
adjusted  to  model  the  stack.  The  membrane  thickness  is  deter¬ 
mined  assuming  that  the  membrane  is  Nation  117.  In  the  MEA 
model,  this  assumption  is  valid,  because  a  membrane  hydration 
model  and  water  transport  mechanisms  that  are  fitted  to  Nation 
117  are  adopted. 

Fig.  5  shows  the  electric  load  current,  while  Fig.  6  shows  the 
relationship  between  the  temperature  and  stack  current  varia¬ 
tions.  The  MEA  temperature  as  the  equal  to  the  stack  temperature 
corresponding  to  the  heat  capacity  of  a  PEFC  stack.  When 
the  load  current  increases  from  50  to  250  A,  as  in  Fig.  5,  the 
stack  temperature  increases  above  75  °C.  At  the  same  time, 
the  convective  heat  transfer  between  the  MEA  and  the  channel 
increases,  and  the  increment  of  convective  heat  transfer  crucially 
affects  the  channel  temperature.  The  anode  and  cathode  chan¬ 
nel  temperatures  rapidly  rise.  When  the  load  current  decreases 
from  250  to  150  A,  as  in  Fig.  5,  at  6000  s,  the  stack  temperature 
rapidly  decreases  to  around  57  °C.  This  decrease  in  temperature 
occurs  because  of  the  reduction  of  heat  generation  due  to  the 


Time  [s] 

Fig.  6.  Temperature  variation  with  anode  inlet  RH  100%  and  cathode  inlet  RH 
100%  at  65  °C. 


Fig.  7.  Total  reaction  energy  rate,  output  electric  power,  and  heat  dissipation 
with  anode  inlet  RH  100%  and  cathode  inlet  RH  100%  at  65  °C. 


reaction  and  to  irreversible  losses  such  as  Ohmic  loss,  activation 
and  concentration  losses. 

The  reduction  of  power  generation  by  the  electrochemical 
reaction  occurs  at  6000  s  in  Fig.  7,  and  the  effects  of  irre¬ 
versible  losses  are  discussed  later.  As  the  temperature  of  the 
stack  changes,  the  temperature  of  the  cooling  water  exiting  the 
stack  also  changes  due  to  difference  between  the  inlet  cooling 
water  and  stack  temperature.  The  convective  heat  transfer  by  the 
cooling  water  accounts  for  a  large  portion  of  power  dissipation 
compared  with  other  factors  such  as  convective  heat  transfer  to 
ambient  air  and  power  converted  to  electricity.  The  convective 
heat  transfer  by  water  cooling  is,  in  particular,  more  significant 
than  the  electric  power  output  from  the  PEFC  at  high  load  cur¬ 
rents.  These  facts  show  that  the  increase  of  irreversible  losses 
reduces  the  efficiency  of  the  PEFC  under  higher  current  loads. 

If  the  load  current  increases,  the  cell  voltage  decreases  at 
3000  s  (Fig.  8).  The  unit  cell  voltage  decreases  after  3000  s. 
Although  activation  losses  become  smaller  due  to  increasing 
temperature,  the  Ohmic  and  concentration  losses  increase  and 
they  have  the  main  effect  on  the  decrease  of  the  cell  voltage  from 
0.73  to  0.44  V.  After  4000  s,  when  the  load  current  decreases, 
the  voltage  output  increases  to  0.59  V  because  of  the  reduc¬ 
tion  in  electric  load  current,  which  reduces  the  Ohmic  losses 
and  concentration  losses.  In  the  Fig.  8,  the  open-circuit  volt¬ 
age  of  the  unit  cell  varies  little  with  changes  in  the  electric 
load  current:  This  change  is  negligible  compared  with  that  of 
other  losses.  The  change  in  Ohmic  loss  during  the  period  of 
constant  load  current  (4000-6000  s)  is  definitely  larger  than  the 
changes  in  activation  and  concentration  losses.  Ohmic  loss  is 
sensitive  to  the  membrane  water  content  which  is,  in  turn,  mainly 
affected  by  the  relative  humidity  of  inlet  gas  and  the  membrane 
temperature. 

Fig.  9  shows  the  transition  of  membrane  water  content  dur¬ 
ing  operation  using  fully  humidified  feed  gas.  The  membrane 
water  content  is  determined  by  the  amount  of  water  surrounding 
the  membrane  and  the  membrane  temperature.  In  the  previous 
section,  the  relationship  between  membrane  water  content  and 
membrane  temperature  was  presented  as  a  function  of  water 
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Fig.  8.  Unit  cell  voltage  and  voltage  losses  with  anode  inlet  RH  100%  and 
cathode  inlet  RH  100%  at  65  °C. 

activity,  which  is  determined  by  the  membrane  temperature  and 
the  water  vapour  pressure  on  both  sides.  When  the  PEFC  oper¬ 
ates  at  a  load  current  of  50  A,  the  membrane  temperature  is 
below  45  °C,  so  that  the  amount  of  water  produced  and  sup¬ 
plied  is  sufficient  to  hydrate  the  membrane.  When  the  PEFC 
operates  at  high  load  current  (250  A),  the  temperature  of  the 
MEA  rises  to  75  °C,  as  shown  in  Fig.  6.  Thus,  the  membrane 
water  content  becomes  around  14.4  (Fig.  9).  During  this  stage, 
the  membrane  water  content  decreases  as  the  MEA  tempera¬ 
ture  increases,  which  affects  the  reference  state  of  the  water, 
namely,  the  water  saturation  pressure.  After  4000  s,  the  load 
current  is  reduced  to  150  A,  as  shown  in  Fig.  5.  The  decrease 
of  water  production  which  occurs  because  of  the  decrease  of 
load  current  causes  an  under- shoot  of  the  membrane  water 
content  for  a  short  time.  After  this  under- shoot,  the  water  con- 


Fig.  10.  Pressure  of  gas  in  anode  channel  with  anode  inlet  RH  100%  and  cathode 
inlet  RH  100%  at  65  °C. 


tent  rebounds  because  the  MEA  temperature  decreases.  In  this 
model,  the  water  content  shows  the  dynamic  behaviour  deter¬ 
mined  by  the  change  in  membrane  temperature.  The  simulation 
shows  that  suitable  water  and  thermal  management  are  needed 
simultaneously. 

Another  important  operating  condition  is  channel  pressure. 
The  partial  pressures  of  the  species  in  the  channel  are  deter¬ 
mined  by  the  mass  continuity  and  the  Ideal  Gas  Eaw.  In  this 
study,  all  of  gas  species  are  modelled  using  the  Ideal  Gas  Eaw. 
Figs.  10  and  11  present  the  pressure  in  the  anode  and  cathode 
channels,  respectively.  The  change  in  the  total  pressure  in  each 
channel  is  unremarkable  because  it  is  assumed  that  the  flow 
rate  of  the  exhaust  gas  is  proportional  to  the  channel  pressure 
and  that  the  inlet  flow  rate  is  constant.  The  channel  pressure 
changes  slightly  because  of  the  transport  of  water  across  the 
MEA.  On  the  other  hand,  the  partial  pressures  of  species  such 
as  water  vapour,  hydrogen,  oxygen  and  nitrogen  exhibit  rel¬ 
atively  large  variations.  If  the  water  vapour  is  saturated,  the 
phase  of  water  changes  from  vapour  to  liquid  based  on  the 
channel  temperature.  Thus,  an  increase  in  channel  tempera- 


Fig.  9.  Membrane  water  content  with  anode  inlet  RH  100%  and  cathode  inlet  Fig.  11.  Pressure  of  gas  in  cathode  channel  with  anode  inlet  RH  100%  and 
RH  100%  at  65  °C.  cathode  inlet  RH  100%  at  65  °C. 
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ture  causes  an  increase  in  vapour  pressure.  Fig.  10  shows  that 
change  in  of  water  vapour  pressure  in  anode  channel  follows 
a  similar  trend  to  the  change  in  the  anode  channel  tempera¬ 
ture. 

Fig.  1 1  gives  the  pressure  in  the  cathode  channel  and  the  par¬ 
tial  pressures  of  the  involved  species,  namely,  oxygen,  nitrogen, 
and  water  vapour.  In  the  cathode  channel,  the  water  vapour  pres¬ 
sure  varies  with  the  cathode  channel  temperature.  Conversely, 
the  partial  pressures  of  the  other  species  such  as  oxygen,  and 
nitrogen  vary  inversely  with  vapour  pressure.  Although  the  par¬ 
tial  pressures  of  hydrogen,  oxygen,  and  nitrogen  in  the  channel 
are  also  affected  by  the  temperature  (per  the  Ideal  Gas  Law), 
the  amounts  of  these  species  decrease  with  increasing  quantities 
of  exhaust  gas  due  to  the  difference  between  the  channel  and 
ambient  pressures. 

In  the  case  of  water  vapour,  liquid  water  is  present  in  the 
channel  (Fig.  12),  that  is,  the  vapour  pressure  is  determined 
only  by  the  relationship  between  saturation  vapour  pressure  and 
temperature  under  this  operating  condition. 

The  net  water  transport  across  the  channel  is  from  the  cath¬ 
ode  to  the  anode  side.  The  upper  graph  of  Fig.  13  plots  the 
direction  and  amount  of  the  net  water  transport  across  the  mem¬ 
brane.  The  lower  graph  in  Fig.  13  shows  that  during  increase 
in  the  load  current,  the  electro-osmotic  drag  suddenly  increases 
up  to  4  x  10_3kgs_1,  and  the  convective  water  flows  due  to 
the  pressure  gradient  increase  simultaneously  in  the  opposite 
direction.  A  negative  sign  indicates  that  the  water  moves  from 
the  cathode  to  the  anode  channel.  Water  transport  from  the 
cathode  to  the  anode  side  assists  the  water  management  of 
the  fuel  cell,  because  it  prevents  the  drying  out  of  the  anode- 
side  membrane  by  electro-osmotic  drag.  After  an  increase  in 
the  electro-osmotic  drag  occurs,  the  magnitudes  of  the  electro- 
osmotic  drag  and  the  convection  by  pressure  gradients  decrease 
to  around  3  x  10-3  kg  s_1.  The  slow  change  of  the  MEA  tem¬ 
perature  affects  the  membrane  water  content,  which  in  turn 
determines  the  electro-osmotic,  water  transport  and  diffusion 
coefficients.  The  diffusion  of  the  water  is  relatively  small  com¬ 
pared  with  the  electro-osmotic  drag  and  convection  by  the 
pressure  gradient. 
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Fig.  13.  Water  transport  across  MEA  with  anode  inlet  RH  100%  and  cathode 
inlet  RH  100%  at  65  °C. 


The  simulation  was  conducted  again  by  changing  the  inlet 
humidity  condition.  The  relative  humidity  of  the  inlet  gas  was 
set  to  75%  at  a  65  °C  dew  point.  Fig.  14  plots  the  dependence  of 
temperature  on  humidity.  Compared  with  Fig.  6,  the  stack  tem¬ 
perature  increases  to  78.4  °C  at  a  load  current  of  250  A,  due  to 
an  increase  in  Ohmic  loss.  The  increase  in  Ohmic  loss  causes  an 
increase  in  the  heat  converted  from  the  produced  power  due  to  the 
electrochemical  reaction  instead  of  electric  power.  Low  humid- 


Fig.  12.  Liquid  water  in  the  anode  channel  and  cathode  channel  with  anode  inlet  Fig.  14.  Temperature  variation  with  anode  inlet  RH  75%  and  cathode  inlet  RH 

RH  100%  and  cathode  inlet  RH  100%  at  65  °C.  75%  at  65  °C. 
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Fig.  15.  Unit  cell  voltage  and  voltage  losses  with  anode  inlet  RH  75%  and 
cathode  inlet  RH  75%  at  65  °C. 

ity  causes  a  decrease  in  the  membrane  conductivity  and  hence  an 
increase  in  Ohmic  loss,  as  illustrated  in  Fig.  15.  The  lower  graph 
shows  the  increase  in  Ohmic  loss  due  to  low  humidity  as  well  as 
the  time- varying  Ohmic  loss  due  to  the  increase  in  stack  tempera¬ 
ture  under  a  current  load  of  250  A  (Fig.  14).  Under  fixed  humidity 
conditions,  the  membrane  humidity  level  is  higher  at  higher 
temperatures. 

The  water  content  reported  in  Fig.  16  is  smaller  than  that 
in  Fig.  9,  because  the  inlet  gas  has  a  lower  humidity.  In  the 
transient  region,  the  graph  does  not  vary  smoothly  near  a  water 
content  of  14  because  the  relationship  between  water  activity 
and  water  content  is  changed  near  that  point  due  to  Eq.  (28). 
The  relationship  between  water  activity  and  membrane  water 
content  above  the  vapour  saturation  level  is  not  correctly  corre- 
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Fig.  16.  Membrane  water  content  with  anode  inlet  RH  75%  and  cathode  inlet 
RH  75%  at  65  °C. 


Fig.  17.  Unit  cell  voltage  and  voltage  losses  with  anode  inlet  RH  75%,  cathode 
inlet  RH  75%  at  65  °C  and  a  cooling  water  temperature  of  35  °C. 


lated  in  the  literature.  In  this  model,  the  Ohmic  loss  and  stack 
temperature  recursively  influence  each  other  through  the  mem¬ 
brane  water  content  and  water  activity.  To  avoid  an  increase 
in  Ohmic  loss  under  low-humidity  conditions,  the  tempera¬ 
ture  of  the  cooling  water  should  be  controlled.  Fig.  17  shows 
reversibilities  of  the  stack  at  a  cooling  water  inlet  temperature  of 
35  °C. 

Cooling  water  at  a  low  temperature  aids  in  the  water  man¬ 
agement  of  the  stack.  A  low  inlet  cooling  water  temperature 
improves  exhausting  heat  and  regulating  stack  temperature. 
At  low  temperatures,  a  relatively  small  amount  of  water  feed 
through  the  inlet  gas  is  sufficient  to  hydrate  the  membrane, 
thus  increasing  conductivity.  The  Ohmic  loss  is  therefore  less 
than  that  when  the  cooling  water  is  at  40  °C  (Fig.  17).  In  order 
to  operate  the  PEFC  stack  under  optimum  conditions,  external 
humidification  and  cooling  water  temperature  should  be  applied 
simultaneously. 

5.  Summary  and  conclusions 

In  this  study,  a  PEFC  model  structure,  which  consists  of  three 
control  volumes,  is  developed.  The  three  control  volumes  are 
the  anode,  the  cathode  and  the  ME  A  control  volumes.  The  mass 
conservation  law  and  energy  conservation  law  are  applied  to 
each  control  volume  in  order  to  model  the  flow  dynamics  and 
thermal  behaviour  of  the  species  therein. 

The  simulation  exhibits  reasonable  behaviour  for  the  mass 
and  energy  flows  within  the  PEFC.  However,  some  simu¬ 
lation  parameters  from  other  literature  are  used  [3-5,8].  In 
future  studies,  these  parameters  will  be  adjusted  using  experi¬ 
mental  data.  This  lumped  parameter  model  is  mathematically 
simple  but  describes  the  essential  states  of  the  system  that 
affect  the  performance  of  the  PEFC,  such  as  temperature, 
pressure  and  humidity.  This  modelling  study  is  therefore  use¬ 
ful  for  understanding  interactions  among  control  volumes 
and  for  designing  model-based  controllers  for  PEFC  sys¬ 
tems. 
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